Order parameter expansion study of synchronous firing induced by quenched noise in 

the active rotator model 
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We use a recently developed order parameter expansion method to study the transition to syn- 
chronous firing occuring in a system of coupled active rotators under the exclusive presence of 
quenched noise. The method predicts correctly the existence of a transition from a rest state to a 
regime of synchronous firing and another transition out of it as the intensity of the quenched noise 
increases and leads to analytical expressions for the critical noise intensities in the large coupling 
regime. It also predicts the order of the transitions for different probability distribution functions 
of the quenched variables. We use numerical simulations and finite size scaling theory to estimate 
the critical exponents of the transitions and found values which are consistent with those reported 
in other scalar systems in the exclusive presence of additive static disorder. 



I. INTRODUCTION 

In some cases, a dynamical system with many variables 
depends on a set of parameters which, although fixed 
in time, are randomly distributed according to a given 
probability distribution. The outcome of the system, al- 
though deterministic, depends on the actual realization 
of the set of parameters. The influence of this so-called, 
depending on the context: quenched noise, static disor- 
der, heterogeneity, variability, diversity, impurities, etc. 
has been the subject of many investigations. In the last 
years, some emphasis has been put in identifying those 
situations in which the presence of the quenched noise 
induces some sort of macroscopic ordering, such as phase 
transitions [1] , patterns 2] ; improves the global response 
to an external forcing [3( or enhances synchrony of firing 
units Q . 

Due to the complexity of the problem, the analytical 
treatments are usually very difficult to be carried out in 
full detail and most results rely on extensive numerical 
simulations. However, a recently introduced technique 
named "order parameter expansion" |5l-ll0f offers a sim- 
ple approximate way of analyzing the effect of the random 
quenched terms in the dynamical equations. The approx- 
imation scheme allows the reduction of the large number 
of coupled differential equations for the microscopic vari- 
ables to just a few effective equations for the relevant 
macroscopic dynamical variables: the mean values and 
dispersions from the mean. 

It is the purpose of this paper to apply this technique 
to the study of the active rotator model |ll| - tl3| under 
the influence of static disorder in the natural frequen- 
cies. Previous work [4. ll4l - tl6| has shown the paradoxical 
result that intermediate levels of disorder at the micro- 
scopic level induce macroscopic order which manifests it- 
self in a coherent firing of the units. Although very so- 
phisticated treatments of this model do exist leading to 
analytical solutions in some particular cases |17H19| (and 
we will refer to them later in the paper) a particular sim- 
ple analysis was developed in [4| , where the authors used 



an expansion of the dynamical equations of the model up 
to first order in the deviation of the quenched variables 
from their mean value and identified a self-consistency 
equation which had to be solved numerically. The or- 
der parameter expansion used in the present article ex- 
pands consistently this analysis up to terms of second 
order, thus reaching a higher accuracy. The resulting 
closed system of three differential equations reproduces 
the ordering abilities of quenched noise in this system 
and, furthermore, predicts a sharp transition back into 
the disordered state where no macroscopic order is ob- 
served. 

The article is organized as follows. First, in section [TT1 
we will define the active rotator model and summarize 
its main properties. Macroscopic observables that de- 
scribe the collective behavior are introduced. Then, in 
section Hnl the approximation method is applied and con- 
clusions about steady states are drawn. In section [IV] 
we present numerical results that support the previous 
findings and use the theory of finite-size scaling to de- 
termine some of the critical exponents characterizing the 
transitions. The paper closes with concluding remarks in 
section W\ 



II. MODEL 

Let us consider a system of globally-coupled active ro- 
tators [ll[ , defined by a set of angular variables <pi (i), i — 
1 , . . . , N which evolve according to the dynamical equa- 
tions: 
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C is the coupling constant. The so-called natural fre- 
quencies uii are quenched noise, i.e. random variables in- 
dependently drawn from a probability distribution g(u>) 
with mean (w). The variance of the distribution, a 2 , is a 
measure of the dispersion of natural frequencies amongst 



the oscillators and measures the degree of intrinsic disor- 
der. We will refer to a as the "diversity" . 

For uncoupled units, C = 0, a value of |w,| > 1 results 
in a rotating behavior for fa{t). The actual period of 

rotation is , = and the direction of rotation depends 
y/wf-i 

on the sign of u>i : clockwise if Wj < and ant i- clockwise 
otherwise. If |w»| < 1, then unit i is in an excitatory 
regime. In this case there are two fixed points (one sta- 
ble and the other unstable) located at the two solutions 
of <p* = arcsin(wi). When a perturbation is such that it 
makes variable fa to cross over the unstable fixed point, 
the subsequent dynamics returns to rest again in the sta- 
ble fixed point through a full turn of fa on the unit circle 
(a "spike" or a "pulse"). This is the typical behavior of 
an excitable system |2C| . 

When the coupling is active, C > 0, the dynamics of 
each unit is influenced by the others which act, effectively, 
as a perturbation. As a result, individual spikes can be 
generated. Those spikes can be independent of each other 
or, alternatively, the units might spike with some degree 
of synchrony. It is of interest to characterize the global 
behavior of the system in order to identify the region of 
parameter space for which synchronized spiking occurs. 
To this end one usually defines a complex variable which 



represents the center of mass of all rotators [21 
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Henceforth, (■ ■ ■) denotes an average over the N units. 
The Kuramoto order parameter p = p(t), where the 
overline denotes an average with respect to time, dif- 
ferentiates between fully synchronized {p = 1, i.e. 
fa(t) — tf>j(t),Vi,j) and desynchronized oscillators {p s=s 
0). When p is close to 1, one still needs to distin- 
guish the rest state where all oscillators are equally 
constant in time from the coherent firing regime where 
the units are oscillating synchronously. Amongst other 
possible measures, one can use the order parameter 
introduced by Shinomoto and Kuramoto [14j as £ = 



p(t)e**(*) — p(t)e l *(*) , which is different from zero only 

in the case of synchronous firing. Alternatively, and this 
is the approach followed in this paper, one can measure 
the average angular speed of the time evolution of the 
global phase Vl / (t). In the rest state, \l/(t) is time inde- 
pendent and the angular speed is zero, whereas in the 
coherent firing regime, ^(t) increases with time and the 
angular speed adopts a non-zero value. 

It has been shown that the system of coupled active 
rotators displays a disorder-induced transition from the 
global rest state to synchronized firing 0, H3, LL5| . Higher 
levels of disorder lead the system again into unsynchro- 
nized firing. The disorder can be originated by the exis- 
tence of diversity amongst the natural frequencies [l6[ (as 
it is the case of interest in this paper) , by the presence of 
noise terms in the dynamical equations, by the existence 
of competitive interactions, heterogeneity in the network 



of connectivities [22[ or any other origin. A general the- 
ory to explain this transition has been developed in j^j, 
while an exact treatment in the case of disorder in the 
natural frequencies has been carried out in [18l , ll9| . In the 
next section we present a simple treatment of this prob- 
lem in terms of the order parameter expansion, which 
allows one to derive equations for the macroscopic vari- 
ables as a perturbative expansion, assuming small fluctu- 
ations. This simple approach is able to predict the main 
features observed in the numerical simulations. Further- 
more, it gives access to an analytic expression for the 
critical noise intensities in the large coupling limit. 



III. METHOD 

A. Derivation of the dynamical equations 

As stated in the introduction, our goal is to use the 
order parameter expansion method to obtain evolution 
equations for the global phase fy(t), defined in Eq.®, 
and its fluctuations, defined as suitable moments of the 
variables Ci(t) — fait) — &(i). We first notice that ac- 
cording to this definition and using Eq. ([2]), a short al- 
gebra leads to (e^W) = e J *W (cosej(t) -Msinej(t)) = 
p(t)e 1 ^^ . Since p(t) has to be a real number we find 
that (sinej(t)) = and p(t) = (cose(i)). As a conse- 
quence we can rewrite Eq. (TTJ) as: 

fait) =Ui- sm.(fa(t)) - C (cosej(i)) sine^i) . (3) 



If we now write Si = Ui — (w) for the deviation of 
the local natural frequency from the mean and take 
then the time derivative of @ one can identify real and 
imaginary parts and find the identity 4?(t) (cos ej(t)) = 
( 4>j (t) cos 6j (t) \ . There we substitute fa (t) by Eq. ^ 

and obtain an equation for ^(t) as a function of (cose.,), 
(cos 2 €j), (sine.,- coscj) and (Sjcosej). If we now expand 
these four terms up to second order around e^ = we are 
left with: 



<j = (u) -fi 2 (£) sin <£(£). 



(4) 



Here we have identified the dynamical variable ^(i) = 
1 — — . We determine its dynamics by writing 

rt 2 (*) = -<Cj(t)ej(t)>. using e z (t) = fa(t) - V(t), and 
replacing fa from Eq. © and "J from Eq. (0|. Expand- 
ing the resulting expression up to second order around 
e; = we arrive at: 



n. 



-W(t) +2 (cos* {t) + C) {1 - Sl 2 (t)) , (5) 



where the third dynamical variable Wit) = (ej(t)5j) al- 
lows us to close the set of equations. It obeys dynamics 
given by W(t) = (ej(t)Sj) and is found in the same way 
as above: 



W = a 1 - (cos *(t) + C) W(t) , 



(6) 



where we made use of the definition (8?) = a 2 . 

The set of equations ^ for the global phase and (JSJB]) 
for the fluctuations, is the result of the order parame- 
ter expansion applied to the oscillator ensemble defined 
by Eqs. ([IJ and is the basis of our subsequent analysis. 
The errors are of the order O ((<5"e™)) , n + m = 3, or 
higher. As a consequence, Eq. (IJ) is more accurate than 
the corresponding equation ^ = (uj) /p — sm^ + O ({$])) 
obtained in [J] . Note that this last equation simply iden- 
tifies p as the threshold for excitability. In our case, the 
full stability analysis is more involved as ^2 (t) is consid- 
ered to be a variable of time. In the next subsection we 
will determine the fixed points of the system (|4j6]) and 
their stability. 



B. Phase diagram 

The fixed points (i&*,Q,%,W*) of the system of equa- 
tions (|4][6|) must satisfy: 



(u) = fi^sintf*, 



n* 2 = l- 



w* = 



„2 



cos <J * + C 



(7) 
(8) 

(9) 



Graphically, the coordinates \&* of the fixed 
points correspond to the intersections of the function 
filO^*) sm C^*) with the horizontal line representing (ui). 
As shown in figure [TJ it turns out that, for fixed (uj) and 
C there exist two limiting values of the diversity a c and 
a' c such that two solutions are found whenever a < a c or 
a > <t' c . A linear stability analysis shows that in this case 
the global phase behaves as an excitable system, corre- 
sponding one of the solutions to an stable and the other 
to an unstable fixed point. If, otherwise, a € (a c , cr' c ) 
there will be no fixed points and the global phase will 
rotate in time, signaling the existence of coherence firing 
in the global system. The linear stability analysis also 
shows that the stable and unstable fixed points, found in 
the low and high noise limits, collide and disappear when 
the maximum of the right-hand-side of Eq. J7J coincides 
with (uj) . This is a so-called SNIC (saddle- node on an in- 
variant cycle) bifurcation [2J]. The steady states in the 
macroscopic equations (HHS]) at high and at low values of 
<7 are caused by different underlying microscopic dynam- 
ics: whereas individual rotators are moving at high levels 
of noise, they are all at rest in the low noise limit. 

The phase diagram identifying regions of synchronized 
global firing can be obtained from the existence of solu- 
tions to equations (JTHHJ) as discussed above. In general, 
this has to be performed numerically, but to an arbitrary 
degree of accuracy. Results for the case that the mean of 
natural frequencies is (ui) = 0.97 are shown in figure [2] 
It can be observed that a minimal coupling intensity is 
needed to introduce a possible state of coherent firing. In 
the large coupling limit, C> 1, it is possible to derive 
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FIG. 1: Graphical analysis of the solutions of the equa- 

tion (u) = (l - 2(cos g! +c) . ) sin(**) = n3(**)sin(**) for 

C = 4.0 and a — 0.5,3.0,7.5 (dashed, dotted, dash dotted). 
The horizontal line marks (cj) = 0.97. Note that the line cor- 
responding to a = 3.0 does not cut the horizontal line, thus 
no stable steady state exists for this value of <r, whereas two 
solutions exist for the other values of a. 



analytical expressions for a c and a' c . Neglecting cos^* 
in the denominator of the right hand side of Eq.®, the 
necessary condition | (uj) \ < f^ leads to 



0-c 



2a/1 
2-v/l 



(10) 
(11) 



In this approximation, the width of the interval (a c , a' c ), 
where the system fires synchronously, grows linear with 
C. This means that an intermediate level of coupling is 
needed to support a synchronized firing state. The de- 
pendence on (uj) of the second transition is rather small 
for (uj) ~ 1. The interval collapses for (ui) = 0. The 
resulting approximate phase diagram for large coupling 
values is marked with dashed lines in figure [2] We con- 
clude that the order parameter expansion correctly iden- 
tifies the diversity induced transitions that occur at the 
critical points a c and a' c . As shown in figure [2J it also 
allows the determination of the value of a c with a rea- 
sonable accuracy, although a' c is grossly overestimated, 
when compared against the numerical simulations (see 
section IIVI) or the analytical treatment of [l9( using a 
Gaussian distribution g(uJi) for the natural frequencies. 
As an example, for (ui) = 0.95 and coupling C — 4, the 
numerical solution of Eqs. (JTHHJ) yields a critical noise in- 
tensity of <j c — 1.269, whereas the approximate solution, 
Eq, (fTUf yields a c = 1.265. This is to be compared with 
the value av = 1.272 obtained from the exact treatment 
given in [19j based on recent developments by Ott and 
Antonsen Il7l.ll8l. 





FIG. 2: Phase diagram: Fixed points of Eqs. (f7ll9|) (for (u>) = 
0.97) exist below and above the two continuous lines (gray 
region) . In between no fixed points exist and the global phase 
rotates, i.e. the individual rotators oscillate in a coherent 
manner. Approximate values of the critical diversity for C S> 
1 according to Eqs. (|10lll[) are plotted with dashed lines. 
Symbols show values taken from numerical simulations of the 
set of Eqs. ([TJ with a Gaussian distribution. The dotted line 
marks the approximate solution of the critical diversity given 
in []]|. 



FIG. 3: Phase diagram for the exponential distribution of 
natural frequencies which satisfies a — (u)}. As in figure [2] 
fixed points of Eqs. (17I9[I exist below and above the two con- 
tinuous lines (gray region). Approximate values of the critical 
diversity for C>1 according to Eqs. (|12II13|) are plotted with 
dashed lines. Symbols show numerical simulations with ex- 
ponential distribution. 



From the microscopic point of view, one could argue 
that the destruction of coherence at high noise values is 
due to the coexistence of individual oscillators rotating 
at opposite directions, as they would certainly be present 
for many general distributions g(w) of natural frequen- 
cies. However, the only requirements we have made on 
the distribution g (uS) is that its first and second moments 
are well defined. Therefore, according to our treatment, 
the existence of elements rotating in both directions can 
not be the only responsible for the transitions. To an- 
alyze this issue, we have considered that the individual 
frequencies were drawn from an exponential distribution 
g(uji) — qt^I^i I (a;), for Wj > such that all natural fre- 
quencies Ui would be positive. In this case the variance 
a 2 and the mean (ui) are not independent of each other, 
as they satisfy a = (ui) and there is only one parame- 
ter in the distribution. Replacing a = (uj) in Eqs. (|T0|) 
and (fTTj) we obtain 



C (V C 2 + 2-CJ 
C ( VC 2 + 2 + c) , 



(12) 
(13) 



as the limits of the zone for which synchronized firing 
exists. The phase diagram for this exponential distri- 
bution has been plotted in figure [3] As it is a special 
case of the general distributions considered above, the 
qualitative image is the same: an intermediate value for 
the intensity of the quenched noise is required to induce 



a state of coherent firing, while a too high intensity de- 
stroys it. As in the case of Gaussian distribution of natu- 
ral frequencies, the qualitative picture agrees with the ex- 
act treatment and the numerical simulations. The lower 
critical point a c is also given with a reasonably degree of 
accuracy, but the upper critical point is overestimated, 
again compared with the numerical simulations or the 
analytical treatment of (l9| . 

To end this section, we note that, for a general distri- 
bution g(u>), a very large diversity satisfying a > a' c will 
never induce another SNIC bifurcation into a new state 
of coherent firing. With this observation one would ex- 
pect that distributions g(oj) with infinite variance, as is 
the case for a Lorentzian distribution, would never show 
a regime of synchronized firing. This is in agreement 
with the detailed theory of [19| only for (w) < 1. In 
the case (u) > 1, however, the complete theory predicts 
that oscillators rotate coherently for low diversity and 
incoherently for high diversity. 

In summary, and in agreement with more involved 
treatments of the coupled active rotator model, the order 
parameter expansion scheme predicts a transition into co- 
herent firing and out of it, induced by the exclusive pres- 
ence of quenched noise. The only assumptions we made 
on the frequency distribution to derive the results are 
the existence of well defined first and second moments. 
In the following section we present numerical simulations 
of the full system, Eqs. ([1]). 



IV. NUMERICAL SIMULATIONS 

In the previous section we demonstrated that for very 
low and very high values of a the system §%HJ§§ is in a 
steady state characterized by "J = (I2 = W = 0, whereas 
for intermediate values the global phase VP is not con- 
stant. This reproduces, in a simple manner, the predic- 
tion of the existence of this intermediate level of disorder 
for which the system fires synchronously and shows the 
validity of the order parameter expansion applied to this 
model. In this section, we will present results of numeri- 
cal simulations of the full system of coupled equations ([1]). 
Our goal is to show that the transitions occurring at a c 
and a' c show all the features of true phase transitions and 
can be characterized, besides by the vanishing of the or- 
der parameter, by a divergence of the fluctuations. The 
divergence, as usual, is smeared out by finite size effects 
and it is possible to carry out an analysis in terms of 
finite size scaling with the number N of rotators [24| . 
Furthermore we want to compare the macroscopic be- 
havior of systems with symmetrically distributed natural 
frequencies and systems with only positive frequencies. 
Namely, Gaussian distributions are used in the first case 
and exponential distributions in the second. 

As order parameter, to, quantifying the collective firing 
regime we have chosen the time average of the slope of the 
global phase to — ^>. This is expected to vanish for small, 
a < <7 C , and large a > a' c disorder and be non-zero in be- 
tween. In the figures we plot the ensemble average ((to)) , 

and the normalized fluctuations y ~ — ^ \((m 2 )) — ((to)) 2 ] , 

where ((• • •)) denotes an ensemble average over realiza- 
tions of the random noise terms and initial conditions. 
We present separately the results for Gaussian and for 
exponentially distributed frequencies. 
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FIG. 4: Simulation results for Gaussian distributed un's: 
a) The order parameter ((m)) for C — 4.0 and various values 
of (uj). The location of the second transition changes little by 
small variations in (uj) if it is close to one. Simulations were 
done with N = 102400. b) Ensemble fluctuations (performed 
over 1000 realizations of the quenched noise variables and ini- 
tial conditions) of the order parameter increase with system 
size «w) = 0.97 and C = 4.0). 



A. Gaussian distributed w-j's 

The natural frequencies Wj are drawn from a Gaussian 
distribution of mean (to) and variance a 2 . In figure 0J, 
we present the results for different values of the mean 
frequency (ui) as a function of the noise intensity a. One 
can see that for small a the order parameter ((to)) van- 
ishes or, equivalently, that the global phase is constant 
indicating that all oscillators are in the rest state. When 
reaching the critical value a c , the global phase ^{t) starts 
to rotate, i.e. ^(t) ~ to ^ 0. This is the regime of syn- 
chronized firing where a macroscopic fraction of the oscil- 
lators fire in synchrony. Increasing the diversity over the 
second critical value a' c , the global phase \& is constant 
again. This is the phase where all units fire in a desyn- 
chronized manner. As predicted by Eq. (|11[) the second 
transition is relatively constant regarding small changes 
in (u>) when it is close to one. 

The precise numerical determination of the location of 
the transition points a c and a' c is hindered by the finite 



size effects. We have found that the location of the maxi- 
mum of the fluctuations of to can give us a good estimator 
of the transition points, as it is relatively constant with 
system size, see figure|4}D. The results for different values 
of the coupling strength C are indicated with symbols in 
the phase diagram, figure [5J The first transition is pre- 
dicted with high accuracy whereas the second transition 
is highly overestimated by the order parameter expan- 
sion. Another feature predicted by the order parameter 
expansion, namely the existence of a minimal coupling 
necessary for inducing coherent firing, is indeed observed 
in the simulations. 

In the vicinity of both transitions at a c or a' c , the en- 
semble fluctuations x of the order parameter diverge with 
system size. As figure EJa) shows, the maximum value 
Xmax(N) scales as N c with c = 0.65 ± 0.03 at the first 
transition and c = 0.61 ± 0.07 at the second. Interest- 
ingly enough, the values of the critical exponent at both 
transitions seem to be consistent with the value c = 2/3 
observed in a phase transition induced by quenched noise 
in a Ginzburg-Landau model |10| . It turns out that the 
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FIG. 5: Finite size analysis for Gaussian distributed frequen- 
cies with (u>) = 0.97 and C — 4: a) Linear fits of maximal 
fluctuations yield log (xmax) ~ clog (N) with c = 0.65 ± 0.07 
for the first transition (circles) and c = 0.6 ±0.1 for the sec- 
ond (squares), b) Rescaled order parameter ((m(<r, N)))N b ' 2 
collapses as a function of eN b with exponent b = 1/3, (N — 
12800, ...,204800). 



full dependence of N and a at both transitions can be 
fitted using standard finite-size-scaling theory [24|, [25[ as 
(m(a,N))) = N-»/ 2 f m (eN b ) and x (tr,N) = N c f x (eN b ) 
with e = 1 — it/iTc or e = 1 — cr/cr^ and being / m and 
f x suitable scaling functions different for each one of the 
transitions. Our numerical results are not sufficiently 
precise to allow an accurate determination of the expo- 
nent b, but reasonable scaling collapse of the data, see 
figure [5jb), is achieved using b — 1/3, as suggested by 
the analogy with the Ginzburg-Landau model mentioned 
before. 



B. Exponentially distributed Ui's 

The probability distribution function for the natural 
frequencies is g(u>i) = e - "'/^/ (ui) for Wj > 0. As men- 
tioned before, this distribution has only one parameter as 
the standard deviation is equal to the mean a = (u) . It is 
chosen such that all rotators have natural frequencies in 
the same, anti-clockwise, direction. As shown in figure [5] 
we find the same dynamical regimes as a function of the 
disorder a than in the case of an arbitrary distribution. 
This is in accordance to the theoretical predictions dis- 
played in figure [3] the transition into coherent firing is 
rather constant and happens around <t c ~ 1, the interval 
grows with rising coupling strength and a minimal C is 
needed to induce coherent firing. Again the second tran- 
sition is overestimated. As before, we use the maximum 
of the fluctuations in the order parameter (see figure [Sb 
for the case of the transition at <r c ) to estimate values 
for the critical noise intensities and annotate them in the 
corresponding phase diagram (figure [3]). 

The first transition, into coherent firing, is marked by 
diverging fluctuations (for a particular case, C — 5.0, see 
figure [Bb) which scale in the same way with system size 
N as we have seen in the Gaussian case (see evidence 
in figure [7]). However, in stark contrast to Gaussian dis- 



a) 




300 



200 



100 



b) 



N= 12800 

. N=25600 

A N=51200 

; \ N= 102400 




0.99 



1.01 



1.03 



FIG. 6: Simulation results for exponentially distributed uii's: 
a) The order parameter is non-zero for finite disorder. Fluctu- 
ations show that the first transition takes place around a ~ 1, 
as Eq. (|f2[) predicts for large C. b) Ensemble fluctuations (for 
C = 5.0) diverge at the first transition for increasing N . 



tributions, the simulations with exponential distributed 
frequencies give strong evidence that the transition into 
asynchronous firing is now of first order. We compare 
the histograms of steady states for 1000 noise realiza- 
tions around both transitions in figure [H At the first 
transition (left column) the distribution broadens at the 
critical disorder and moves continuously to higher values. 
On the contrary the equilibria near the second transition 
are narrowly distributed around zero, or around the non- 
zero value in the ordered state, typical for a first order 
transition (right column). 

It turns out that the order parameter expansion de- 
veloped in the previous section predicts that the second 
transition into asynchronous firing occuring at a = a' c . 
is of second order for the Gaussian distribution and of 
first order for the exponential distribution of frequencies. 
The results of the numerical integration of the system 
of equations (OHJ) for selected sets of parameters ((w), 

<7, C) for the mean phase velocity \1/ are plotted in fig- 
ures O (Gaussian) and [TU] (exponential). It is evident the 

jump of vfr at the second critical point a' c in the case of 
the exponential distribution whereas it is continuous for 
the Gaussian distribution. The first transition to syn- 
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FIG. 7: Finite size analysis for exponentially distributed fre- 
quencies (C = 5): a) Linear fit of maximal fluctuations yields 
c = 0.65 ± 0.02. b) Rescaled order parameter with scaling 
exponent b — 1/3. 



chronized firing at a — a c is predicted to be continuous 
independently of the distribution of frequencies. 



V. CONCLUSIONS 

We have used the order parameter expansion to ap- 
proximate the dynamics of the global phase in systems 
of coupled active rotators under the influence of quenched 
disorder. The method leads straightforwardly to a sys- 
tem of three differential equations easier treatable than 
the full system and more accurate than other approxima- 
tions used in previous works. In agreement with exact re- 
sults for the full system, the global phase of the reduced 
system can undergo a transition from a rest state into 
a rotating regime and back into a rest state, when sub- 
jected to increasing diversity. In the rest states, ^(i) is 
time independent and the angular speed is zero, whereas 
in the intermediate regime of coherent firing, \I/(t) in- 
creases with time and the angular speed adopts a non- 
zero value. Our treatment allows us to give analytic ex- 
pressions for the critical disorder values in the limit of 
large coupling. We have seen that the first transition is 
predicted to a high degree of accuracy whereas the second 
is highly overestimated. 

We have used numerical simulations to show that the 
ensemble fluctuations of the order parameter diverge at 
the transition points. The simulations with Gaussian dis- 
tributed frequencies show continuous transitions, both in 
and out of the coherent firing state, but if the frequen- 
cies are distributed according to an exponential distri- 
bution (and therefore the mean and variance are varied 
simultaneously) then the destruction of the ordered state 



is achieved through a first order transition. The order 
parameter expansion scheme predicts this distinction of 
the transitions. A finite-size scaling analysis of the nu- 
merical simulations data indicate that the critical expo- 
nents of the transitions are consistent with those found 
in the athermal Ginzburg-Landau model with additive 
quenched noise. 

First transition Second transition 

o=0.9 o=4.7 

1000--- — i — i — i— i 1000 r 
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FIG. 8: Histogram of 1000 steady states from simulations of 
equations fTJ with exponentially distributed frequencies. At 
the transition into coherent firing (left column) the values 
are distributed around one single value, broadening near the 
critical disorder. The destruction of the ordered state is a 
first order transition (right column), the values are distributed 
narrowly around zero or the non-zero value, C = 5.0. 
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